Variational Monte-Carlo calculation of the nematic state of the two-dimensional 

electron gas in a magnetic field 
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We use a Jastrow-Slater wave function with an elliptical Fermi sea to describe the nematic state 
of the two-dimensional electron gas in a magnetic field and the Monte Carlo method to calculate a 
variational energy upper bound. These energy upper bounds are compared with other upper bounds 
describing stripe-ordered ground states which are obtained from optimized Hartree-Fock calculations 
and with those which correspond to an isotropic ground state. Our findings support the conclusions 
drawn in our previous study where the Fermi-hypernetted chain approximation was used instead of 
the Monte Carlo method. Namely, the nematic state becomes energetically favorable relative to the 
stripe-ordered Wigner crystal phase for the second excited Landau level and below a critical value 
of the layer "thickness" parameter which is very close to its value in the actual materials. 

PACS numbers; 73.43.Cd,73.43.Lp 
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I. INTRODUCTION 

The measurements of Lilly et al^i and Du et ali^ 
reveal strong anisotropic transport properties of the 
two-dimensional electron gas (2DEG) for the half-filled 
Landau-level system for high Landau-levels (LL) and at 
very low temperature. The anisotropic behavior in the 
transport properties is consistent with stripe and bubble 
charge-density-wave phases which were predicted early 
on m Refs. iii by means of Hartree-Fock calculations 
of the 2DEG and were confirmed more recently by numer- 
ical studies of systems with up to 12 electrons^il. How- 
ever, Fradkin et al4 have challenged this interpretation 
and suggested that the anisotropic transport might be 
due to a possible nematic phase of the 2DEG in a mag- 
netic field. This idea finds support in the good compar- 
ison between the results of the temperature dependence 
of the anisotropy of the resistivity obtained by means 
of a Monte Carlo simulation of the nematic phase^ with 
that which has been experimentally observed. In addi- 
tion, the idea is supported by the experiments of Cooper 
et al.^" where an in-plane magnetic field was applied in 
the 2DEG and the results of the experiment were inter- 
preted on the basis of the presence of a nematic state; 
further support of the idea is provided by the fact that 
the theoretically estimated transition temperature from 
an isotropic to nematic phase^^ is of similar magnitude 
as the experimentally determined temperature at which 
the on-set of the anisotropic transport occurs. 

Rather recently we have presentedi^ a variational cal- 
culation of the nematic state as ground state of the 
half-filled Landau-level system in a magnetic field based 
on an ansatz ground state wavefunction proposed by 
Oganesyan et al.^'^ which is of the Jastrow-Slater form 
and is given by the following expression: 
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where Pq is the projection operator onto the lowest LL, 
(t>ki^i) are two-dimensional (2D) plane-wave states. Here, 
Zj = Xj + iyj is the complex 2D coordinate of the j elec- 
tron. This wavefunction is a Jastrow correlated Slater 
determinant with Jastrow part similar to the Laughlin 
stated. This ground-state wavefunction has the same 
form as the form proposed by Rezayi and Read^^, how- 
ever, the single-particle momenta form an elliptical Fermi 
sea as opposed to the circular Fermi sea. There is a 
broken-symmetry parameter which is the ratio a = ki jhi 
of the semi-major k\ and semi-minor ki axes of the ellip- 
tic Fermi sea. Using this wave function to describe the ne- 
matic state we had carried out a variational study of the 
half-filled system using the so-called Fermi-hypernetted- 
chain (FHNC) approximation^^. 

The results of the above mentioned variational calcula- 
tion indicate that there is a certain value of the parameter 
A (A is proportional to the 2DEG layer thicknessi^) be- 
low which the nematic state is energetically favorable as 
compared to the isotropic and the stripe-ordered ground 
states for the second excited LL. It is interesting to note 
that this critical value of A is very close to the value of A 
which can be estimated based on the actual experimental 
conditions which are applicable for the case of the data 
by Lilly et alii and by Du et ali^. However, one of the 
weak points of the above described variational study is 
the fact that the FHNC approximation is plagued by an 
unknown-size error and the results cannot be improved 
in a controlled manner-^,. Therefore, there is a need to 
check the validity of these results and conclusions using 
the variational Monte Carlo method and this task is un- 
dertaken in the present work. 

There is a different variational approach to the prob- 
lem of a broken rotational state of the half-filled LL 
introduced by Ciftja and Wexler— . They have used 
the Fermi-hypernetted-chain (FHNC) approximation to 
study a broken rotational state of the half-filled LL where 
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the symmetry-breaking parameter was introduced in the 
correlation part of the wavefunction as {zi — Zj)'^ — > 
(zi — Zj — a){zi — Zj + a), and they used the standard 
single-particle determinant with a circular Fermi sea. 

In this and in the work of Ref. we considered the 
unprojected wavefunction of the nematic state. The ad- 
vantage of this simplified version is that it has a Jastrow 
form with a Slater determinant so it can be applied di- 
rectly with FHNC and it allows us to study large-size 
systems using the variational Monte Carlo method. The 
paper is organized as follows: In the following Section 
we discuss the formulation and the procedure; in Sec. IIIII 
we present the results and we compare them with those 
obtained for the case of a stripe-ordered state and the 
isotropic state. In Sec. IIVI we summarize the conclusions 
of the present calculation. 



II. METHOD 



For a value of N to be acceptable the number of states, 
i.e., the number of pairs {nx,ny) satisfying the above 
inequality should be equal to A''. For example, for a = 1, 
N can be 1, 5, 9, 13, 21, 25, 29, 37, 45, . . .; for a = 2, they 
can be 1,3,7,11,15,17,21,... 
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We have adopted the toroidal geometry of a square 
with periodic boundary conditions. This geometry has 
the advantage of naturally adapting to the nematic and 
isotropic state wavefunction. There are several steps in 
applying the MC approach for this problem. First, as 
part of the wavefunction of nematic state we construct a 
Slater determinant of plane waves characterized by mo- 
mentum vectors which lie inside an elliptical Fermi sea. 
Second, since the pseudo-potential is ln{r), which is a 
long-range interaction, we need to to take into account 
all periodic image charge interactions. One of the meth- 
ods to do this is the Ewald summation technique. In 
subsection IV Al of the appendix we describe the Ewald 
summation technique for the case of toroidal boundary 
conditions and the ln{r) interaction. In the present sec- 
tion we will discuss our implementation of the MC to 
study the nematic state. 

Given a value of a there are definite values of the num- 
ber of particles N which correspond to a closed shell. 
These definite values of N are calculated as follows. The 



occupied states characterized by kx 
following condition: 



ky must satisfy the 



^ + <1, 



(2) 



where fci and ^2 are the major and minor axis of the 
Fermi sea and given by: 



ki 



k. 



(3) 
(4) 



where p is the uniform particle density of the system. For 
a finite system of size L x L, kx = UajA/c and ky = UyAk 
where Afc = 2tt/L and nx,ny G Z. So one can deduce 



the conditions for r 
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FIG. 1: Occupied states for the nematic state with a = 2 
(top) and a = A (bottom) for the case of 89 particles. 



In Figs, m we present two examples of closed shell 
which correspond to a = 2 and 4 for = 89. Notice that 
with anisotropy parameter a = kx/ky > 1^ the occupied 
states (i.e those satisfying equation [5]) are anisotropically 
distributed along the preferred kx axis. In our MC calcu- 
lation wc will use these cases as well as larger size systems 
up to 145 particles. 

We follow the Metropolis MC scheme for sampling the 
wavefunction where the ratio needed between the new 
and the old wavefunction is: 



'4'{foid) 



exp{u{f 



u{fold)) 



Det{e' 







L»ei(e*'=-''°") 



(6) 



where u(f) is the periodic pseudo-potential which is de- 
rived in the subsection (|V Ap of the appendix. To carry 
out the calculation of the ratio between the Slater deter- 
minant of the new configuration and the old configuration 
efficiently, we use the inverse updating technique devel- 
oped by Ceperley et alf^^. We found that the number 
of MC steps needed for "thermalization" is of the order 
of 10^ and we use of the order of 2 x 10^ MC steps to 
calculate averages of the distribution function. 
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The potential energy of the high LL can be expressecU^ 
. .a the pair distribution function of the LLL using the 
single mode approximation discussed in Ref. [13, namely, 



via 



.9^ 



{r)(fr 



(7) 



where the effective potential Vjj:j{r) for Landau level L 
is the convolution of the effective interactiorti^ 



V{r) = e^eVr^ + 



(8) 



with the L-order Laguerre polynomial; namely, it is the 
Fourier transform of: 



eq 



1 2 



(9) 



In the above formula, A is a length scale which charac- 
terizes the confinement of the electron wave function in 
the direction perpendicular to the heterojunctioni^. 

We use the single mode approximation to calculate the 
interaction energy at high LL (equation [7]) and we are 
only interested in obtaining the pair distribution function 
g{r). The kinetic energy advantage of the isotropic phase 
over the nematic phase is calculated in the subsection 
(|VB|) pf the appendix. The approach can be divided 
into the following steps: 

• The pair distribution function for the LLL for dif- 
ferent anisotropic parameters a is calculated. 

• The single mode approximation^" is used to calcu- 
late the interaction energies at a high LL. 

• The kinetic energy for different anisotropic param- 
eters is evaluated (see appendix IV B|) . 

• We compare total energies of the isotropic and ne- 
matic state to determine at what LL the nematic 
becomes energetically favorable. 

• The optimum value of a is determined by minimiz- 
ing the total energy. 

• The HF results which have been reported so 
far— >^>^ correspond to the case of A = 0. There- 
fore, we needed to carry out Hartree-Fock calcu- 
lations following Refs I2lll22l23l for the case of the 
interaction given by Eq. [5] for A ^ 0. The opti- 
mum total energies of the nematic states will be 
compared with those of the stripe states at differ- 
ent values of A for the 2"*^ excited LL to determine 
a critical value of A below which the nematic state 
maybe energetically favorable. 

• The above mentioned critical value of A is compared 
with the value which corresponds to those samples 
used in the experiment^. 

• A comparison of the MC results to the ones ob- 
tained by FHNC will also be presented. 



III. RESULTS 

The pair distribution function g{r) obtained using MC 
integration has important differences when compared to 
g{r) obtained by FHNGi^ as illustrated in Fig. El Thus, 
it is important to obtain the energies of the nematic state 
at high LL by MC and to compare them with those ob- 
tained by FHNC. 

We first compare the interaction energies obtained for 
different values of a > 1 with the potential energy of the 
isotropic state {a — 1) (Figs. OandS]). Notice from Figs. 
[3] and [4] that the potential energy of the isotropic state 
is lower than the potential energy of the nematic state 
for the f** excited LL and LLL for all values of the pa- 
rameter A. The potential energy is calculated with the 
pseudo-potential obtained using the Ewald sum as dis- 
cussed in subsection (jV A|) of the appendix. Essentially 
the same result is also found with a pseudo-potential ob- 
tained using the Lekner summation technique^. Fur- 
thermore, as shown in appendix IV Bl the kinetic energy 
of the isotropic state is below that of the nematic state, 
and, thus, the total energy of the isotropic state is al- 
ways lower than that of the nematic state. Hence, our 
MC calculation shows that the isotropic state is energet- 
ically favorable as compared to the nematic state for the 
LLL and the 1^* excited LL for all values of the param- 
eter A. Also note that the same conclusion was reached 
using the FHNC techniquc^^ with the same wavefunc- 
tion. These findings solidify the conclusion that for the 
LLL and the 1** excited LL, the isotropic state is more 
stable than the nematic state, which is also in agreement 
with the experimental findings of Refs. [H and U 




FIG. 2: Comparison of the pair distribution function obtained 
by FHNC and MC. 



For the 2"'' excited LL, however, the situation changes 
as illustrated in Fig. [5] The conclusion which can be 
drawn from the comparison of Fig. [5] is that the interac- 
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FIG. 3: Comparison of the potential energy of the nematic 
state calculated for various values of a 7^ 1 as function of A 
with the isotropic state (a — 1) for LLL. 



FIG. 5: Comparison between the potential energy of the ne- 
matic state calculated for various values of the anisotropic 
parameter a 7^ 1 as a function of A and the potential energy 
of the isotropic state {a = 1) at the 2"'' excited LL 




FIG. 4: Comparison of the potential energy of the nematic 



state calculated for various values of a / 1 as function of A pjQ g. Comparison of total energy of the nematic state cal- 
with the isotropic state (a = 1) for the 1" excited LL. culated for various values of the anisotropic parameter a ^ 1 

as functions of A with the isotropic state (q = 1) at the 2"'* 

excited LL 



tion energy of the nematic state is lower than that of the 
isotropic state for all values of A. However, we need to 
compare the total energy of the nematic state with that 
of the isotropic state for the 2"'' excited LL (Fig. [6]). 
From Fig. \6\ we conclude that the nematic state is en- 
ergetically favorable as compared to the isotropic state 
for the 2"'' excited LL for the range of the parameter 
A < 0.4. Note that using FHNC we foundi^ that for 
A < 0.6 the total energy of the nematic state is lower 
than the energy of the isotropic state. In summary, both 
FHNC and MC yield similar conclusions about the sta- 



bility of the nematic state against the isotropic state for 
the 2"'' excited LL. 

In Refs. the stripe-ordered phase was predicted 

based on HF calculations and this ordering can also ex- 
plain the anisotropy observed in the transport properties 
of the 2DEG at low temperature. Therefore, we need 
to investigate the stability of the nematic state against 
the stripe-ordered state as follows. First, we find the op- 
timum energies of the nematic state with respect to the 
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anisotropic parameter a for various values of A. Next, we 
compare these with the optimum energies of the stripe 
state obtained by the HF approximatio n^^'^^i^^ . Calcu- 
lations for the case where A = have been carried out 
in Refs. [2lll2^ and [H. For making a comparison with 
the optimum energy of the nematic state at various val- 
ues of A, we carried out detailed HF calculations for the 
case where A 7^ 0. For the stripe-ordered state, the opti- 
mum energy is obtained by minimizing the energy with 
respect to the uniaxial anisotropy parameter e defined 
in Ref. Fig. [7] shows the comparison of the opti- 
mal energies obtained by MC for the nematic state with 
the optimum (with respect to e) energy for the stripe 
state obtained by HF. Note that, for A > 0.5, the opti- 
mum nematic state is obtained for a = 1, i.e., it is the 
isotropic state. Furthermore, Fig. [7] demonstrates that 
the nematic state is energetically lower than the stripe 
state for the values of A < Ac = 0.37. As discussed 
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FIG. 7: Comparison of the optimal nematic state calculated 
by MC using the pseudo-potential using the Ewald sum with 
the stripe state calculated by HF as function of A. 



earlier the pseudo-potential can be obtained by using ei- 
ther the Ewald or the Lekner summation technique^"*. 
We have also carried out the same calculation using the 
Lekner summation technique and the results obtained are 
in good agreement with those obtained using the Ewald 
summation method. Thus, we can conclude that with 
MC calculation, for A < Ac = 0.37, the energy of the 
nematic state is lower than the stripe-ordered state. 



IV. CONCLUSIONS 

In Fig.[8]the results for the optimum total energy of the 
nematic state obtained with the variational MC method 
is compared with that obtained by FHNC in Ref. [12| and 
with the optimum energy of the stripe-ordered state. The 



critical value of Ac we found from EHNC""^^ is 0.4 which 
is close to the value of 0.37 obtained above by MC. The 
critical value of A corresponding to the sample used in 
experiment^ which was calculated in Ref. |l2|, using the 
conditions of the experiment and sample characteristics, 
is approximately 0.34, which can be below the critical 
value found above. Thus, both MC and FHNC calcula- 
tions indicate that the nematic state might be the state 
observed experimentally for the 2DEG at the heterojunc- 
tion in the samples used in experiment described in Ref[l|. 
There is still a remaining question about the validity 
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FIG. 8: Comparison of the optimum nematic state obtained 
from FHNC and MC with the stripe state obtained from HF 



of our approximation to neglect the projection opera- 
tor in the wavefunction ([1]). However, in both FHNC 
treatments of the probleroiSiii^, where, in addition to ne- 
glecting the projection operator for arguments presented 
there, there was a second rather annoying question (and 
rather straightforward to answer) of the validity of the 
FHNC approximation in evaluating the energy expecta- 
tion value. In the present paper the latter question is 
answered by employing the Monte Carlo method. There- 
fore, we conclude that the present calculations eliminates 
the suspicion that the conclusions drawn in Ref. [l^ might 
be due to an artifact of the FHNC approximation em- 
ployed in Ref. [H. 



V. APPENDIX 

A. Ewald summation technique for the logarithmic 
potential 

The long-range nature of the pseudo-potential ln(r) 
which appears in the exponent of the wavefunction of 
the Jastrow factor in the case of periodic boundary con- 
ditions requires a summation over all periodic image 
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charges. Specifically, the "charge" distribution required 
to give rise to a logarithmic interaction is given as: 



P(^) ^^^i'r-R)+ Pbackground 



(10) 



The two-dimensional (2D) Poisson equation is given by: 

V2$(f) = -27rp(f) (11) 

and its solution in 2D is the logarithmic interaction. We 
need to solve the above equation for a periodic square 
L X L. The idea of Ewald summation is to add around 
each charge an opposite Gaussian charge distribution of 
an appropriately chosen width ^, and, in addition, to 
subtract the same Gaussian charge distribution. Let us 
split p into long-range and short-range portions in the 
following manner: 



p{r) = piir) + P2ir) 



— 



Pbackground 



P2(r) 



E 

R 



5{r-R) 



-ay 

7^ 



(12) 



(13) 



(14) 



01, which corresponds to pi is a short-range potential, 
and, thus, we can calculate (jfi in real space since it con- 
verges very quickly. The other combined charge configu- 
ration, i.e., p2, consisting of the Gaussian and the back- 
ground charge and the corresponding potential is denoted 
by (j)2. Since 02 is a long-range potential it will be cal- 
culated in Fourier space. The solution to each of the 
Poisson's equations for the two charge distributions and 
the corresponding potential is straightforward. We note 
that for our case the "charge" of the particle is = 2m. 
We find 



Mr) 



02 (r) 



4m7r e 

fc#o 



A 



-TO Ei 

R 



P' 



(15) 



(16) 



where k = 2-K/Ln with n E and Ei{t) is the Ex- 
ponential integral function and is defined by: Ei{t) = 

For the Ewald summation, the convergence of (jlSp 
and (flB)) is achieved choosing the width of the Gaus- 
sian charge distribution fj, — 1, the number of cells for 
the sum in pB]) to be 10 and by carrying out the sum in 
momentum space in ([15]) over 200 k-states. 

In order to check the validity of this approach for the 
case of our use of toroidal boundary conditions we cal- 
culated the distribution function and the energy for the 
1/3 {m — 3) case using the expressions (115116^ and our 
results for the energy and distribution function are iden- 
tical to the results of Morf and Halperin^^ who used the 
disk geometry. 

B. Evaluation of kinetic energy of the nematic 
state 



isotropic state. In the single-LL approximation, the ki- 
netic energy is quenched. In addition, the same is true in 
the HF treatment of the stripe, namely, there is no kinetic 
energy due to any correlation factors or operators. While 
this approximation gives a significant difference between 
the potential energy of the isotropic state and the ne- 
matic state, it gives no difference between their kinetic 
energy which is unacceptable because of the difference in 
the geometry of the Fermi sea. We want to estimate this 
difference. We can start with: 



= (V - A)'^F<P + 2 



(V - A)F 



(V - AfF<i> (17) 
V$ + FV2$ (18) 



The first term in the above equation yields: 



huJr 



(V - AyF<P ^ ^i^*, 



(19) 



which is common for all states under our consideration 
so for simplicity we can drop it. The last term is: 



2 ri^P 



So the contribution of the last term is: 



^ 2to* 



(20) 



(21) 



keFFS 



where k € FFS stands for a summation over all vectors 
k in the corresponding filled Fermi sea. The summation 
over the circular Fermi sea in the isotropic case is given 
by: 



1 ^ n^fc2 



h^kl 



N ^ 2m* Am* 

k 



(22) 



and the in the case of the elliptic Fermi sea in the 
anisotropic case the summation is given by 



-E 



2l2 



h^k 

2to* 



r? hi 



4to* 



(23) 



Using the facts that k\ = k\ ■ ki and A:i/fc2 = ck, the 
kinetic energy difference between the isotropic state and 
the nematic state is given as follows: 



l\{KE) = - 



r?k\ (1 - oif 



4to* 



2a 



(24) 
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In this subsection of the appendix, we compute the 
kinetic energy difference between the nematic and the 
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